Calculate the RMS between two observatories



[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/image/work/IMAGE_ARRAY.xlsx');
master_b(1,:) = [];

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);
S = dir('/Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/*.mat');

for i = 1:length(S)-1,
    
    eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(i ,1)) '.mat' ' x_data y_data z_data']);
    
    x_data_1 = x_data;
    y_data_1 = y_data;
    z_data_1 = z_data;
    
    
    for j = i+1:length(S),
        
        metadatacellarray = master_b(i ,:);
        lat = master_a(i ,1);
        long = master_a(i ,2);
        alt = master_a(i ,3);
        eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(j ,1)) '.mat' ' x_data y_data z_data']);
        
        
    end;
end;


% steps

% read x,y, z data
% calculate F
% read x,y,z secular removed data
% find scalar F through IGRF
% compare the above two

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);

nc_fname = '/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/ALE_1995_2012_Absolute.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);

% fit spline to F and remove the secular variations

f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);


L = isnan(f_data);

data_array = f_data(~L);

time_array = time_array_min(~L);


b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,
    
    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    
    v=ppval(time_array_min,sp);
    
else,%data length <=1 year
    
    sp=robustfit(time_array,data_array);
    
    v= sp(1) + sp(2) * time_array_min;
    
end;


f_data_secular_removed = f_data - v';


% now read the secular removed data and calculate F using IGRF

nc_fname = '/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/ALE_1995_2012_Secular_Removed.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);


icount = 1;
for decyear = 1995:0.1:2013,
    [mag_igrf(icount,:)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
    icount = icount + 1;
end;

L = isnan(x_data) | isnan(z_data) | isnan(y_data);

% get values at all points through  interpolation
MainB = interp1(datenum(1995:0.1:2013,1,1), mag_igrf, time_array_min, 'linear');

% get a unit vector in the direction of Main field

Unit_Vector = [MainB(~L,:)./repmat(sqrt(MainB(~L,1).^2+MainB(~L,2).^2+MainB(~L,3).^2),[1,3])];

% find dot product

scalar_f = dot([x_data(~L) y_data(~L) z_data(~L)]',Unit_Vector');


% May 16,2013. Both scalar_f and f_secular_removed are the same



%% load all the secular removed data , find F and save the data.
% While the secular removed f can be calculated from absolute X,Y,Z and
% then spline fitting, this is better (as I found with the above
% experiment)

S = dir('/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/*.nc');
S1 = dir('/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/*.nc');

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);


for i = 1:length(S),
    
    
    
    
    nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/' S(i).name] ;
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
    
    
    icount = 1;
    for decyear = 1995:0.1:2013,
        [mag_igrf(icount,:)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
        icount = icount + 1;
    end;
    
    L = isnan(x_data) | isnan(z_data) | isnan(y_data);
    
    % get values at all points through  interpolation
    MainB = interp1(datenum(1995:0.1:2013,1,1), mag_igrf, time_array_min, 'linear');
    
    % get a unit vector in the direction of Main field
    
    Unit_Vector = [MainB(~L,:)./repmat(sqrt(MainB(~L,1).^2+MainB(~L,2).^2+MainB(~L,3).^2),[1,3])];
    
    % find dot product
    
    scalar_f = nan(size(x_data));
    
    scalar_f(~L) = dot([x_data(~L) y_data(~L) z_data(~L)]',Unit_Vector');
    
    % Now if available calculate the F in the traditional way
    
    nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/' S(i).name(1:14) 'Absolute.nc'];

        
    if (~isempty(dir(nc_fname))),
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
    
    % fit spline to F and remove the secular variations
    
    f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);
    
    
    L = isnan(f_data);
    
    data_array = f_data(~L);
    
    
    %% Now select observatories in pair and compute
    
    time_array = time_array_min(~L);
    
    
    b1 = min(time_array):365:max(time_array)+10;
    
    if length(b1) > 1,
        
        sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
        
        v=ppval(time_array_min,sp);
        
    else,%data length <=1 year
        
        sp=robustfit(time_array,data_array);
        
        v= sp(1) + sp(2) * time_array_min;
        
    end;
    
    
    f_data_secular_removed = f_data - v';
    

    
    % plot the data
    plot(time_array_min, f_data_secular_removed);
    hold on;
    
    end;
    
    plot(time_array_min, scalar_f,'r');
    %save the figure
    
    saveas(gcf,['/Users/manojnair/projects/obs_mag_data/Canada/plots/' S(i).name(1:3) '_F_compare'],'png');
    close all;
    clear mag_igrf;
    
    % save the data
    
    save(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' S(i).name(1:3) '_F'], 'scalar_f');
    
    fprintf('Completed %s \n', S(i).name);
end;










